

% Define useful functions
va = @(vv,zzeta,ell) ( Params.b + vv ) ...
            .* (   1 - Params.barg ...
                 + Params.barg * Params.zl0 / Params.rho ...
                 * Params.kappa * zzeta ./ ( zzeta - 1 ) / ( Params.kappa - 1) ...
               ) ...
            .* ell ;
                                    
% Compute corporate tax credit
% Called TauVA in code but corresponds to corporate tax credit
LogTauVA = zeros(Grids.Nx,1) ;
LogTauVA(Grids.Nx)  = 0 ; % Normalize by no tax/subsidy in the best xation

coef = Params.alpha / ( 1 - Params.alpha ) ; 

i=Grids.Nx ;
while i > 1
      
    dzeta_i = zetaP(i) - zetaP(i-1) ;
    if i < Grids.Nx
        dzeta_i = zetaP(i+1) - zetaP(i) ;
    end
    
    dx_i = (   SbarPrimeOverSbar(zetaP(i)) ...
             + log( Params.B0 / ( Params.b + vP(i) ) ) ...
           ) * dzeta_i ;
    
    LogTauVA(i-1) = LogTauVA(i) - coef * ( - dx_i ) ; 
    
    i = i - 1 ;

end

tauVA = exp( LogTauVA ) ;
